Pade- Type Model Reduction of Second-Order 
and Higher- Order Linear Dynamical Systems 



Roland W. Freund 

Department of Mathematics, University of California at Davis, One Shields 
Avenue, Davis, CA 95616, U.S.A. 
f reund@math . ucdavis . edu 

Summary. A standard approach to reduced-order modeling of higher-order linear 
dynamical systems is to rewrite the system as an equivalent first-order system and 
then employ Krylov-subspace techniques for reduced-order modeling of first-order 
systems. While this approach results in reduced-order models that are characterized 
as Pade-type or even true Pade approximants of the system's transfer function, in 
general, these models do not preserve the form of the original higher-order system. 
In this paper, we present a new approach to reduced-order modeling of higher-order 
systems based on projections onto suitably partitioned Krylov basis matrices that 
are obtained by applying Krylov-subspace techniques to an equivalent first-order 
system. We show that the resulting reduced-order models preserve the form of the 
original higher-order system. While the resulting reduced-order models are no longer 
optimal in the Pade sense, we show that they still satisfy a Pade-type approximation 
property. We also introduce the notion of Hermitian higher-order linear dynamical 
systems, and we establish an enhanced Pade-type approximation property in the 
Hermitian case. 



1 Introduction 

The problem of model reduction is to replace a given mathematical model of 
a system or process by a model that is much smaller than the original model, 
yet still describes — at least approximately — certain aspects of the system or 
process. Model reduction involves a number of interesting issues. First and 
foremost is the issue of selecting appropriate approximation schemes that 
allow the definition of suitable reduced-order models. In addition, it is often 
important that the reduced-order model preserves certain crucial properties 
of the original system, such as stability or passivity. Other issues include 
the characterization of the quality of the models, the extraction of the data 
from the original model that is needed to actually generate the reduced-order 
models, and the efficient and numerically stable computation of the models. 

In recent years, there has been a lot of interest in model-reduction 
techniques based on Krylov subspaces; see, for example, the survey pa- 
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pers jSJ 01 IHj ■ The development of these methods was motivated mainly 
by the need for efficient reduction techniques in VLSI circuit simulation. An 
important problem in that application area is the reduction of very large-scale 
RCL subcircuits that arise in the modeling of the chip's wiring, the so-called 
interconnect. In fact, many of the Krylov-subspace reduction techniques that 
have been proposed in recent years are tailored to RCL subcircuits. 

Krylov-subspace techniques can be applied directly only to first-order lin- 
ear dynamical systems. However, there are important applications that lead 
to second-order, or even general higher-order, linear dynamical systems. For 
example, RCL subcircuits are actually second-order linear dynamical systems. 
The standard approach to employing Krylov-subspace techniques to the di- 
mension reduction of a second-order or higher-order system is to first rewrite 
the system as an equivalent first-order system and then apply Krylov-subspace 
techniques for reduced-order modeling of first-order systems. While this ap- 
proach results in reduced-order models that are characterized as Pade-type 
or even true Pade approximants of the system's transfer function, in general, 
these models do not preserve the form of the original higher-order system. 

In this paper, we describe an approach to reduced-order modeling of 
higher-order systems based on projections onto suitably partitioned Krylov 
basis matrices that are obtained by applying Krylov-subspace techniques to 
an equivalent first-order system. We show that the resulting reduced-order 
models preserve the form of the original higher-order system. While the re- 
sulting reduced-order models are no longer optimal in the Pade sense, we 
show that they still satisfy a Pade-type approximation property. We further 
establish an enhanced Pade-type approximation property in the special case 
of Hcrmitian higher-order linear dynamical systems. 

The remainder of the paper is organized as follows. In Section [2] we re- 
view the formulations of general RCL circuits as first-order and second-order 
linear dynamical systems. In Section |31 we present our general framework for 
special second-order and higher-oder linear dynamical systems. In Section QJ 
we consider the standard reformulation of higher-order systems as equivalent 
first-order systems. In Section |S1 we discuss some general concepts of dimen- 
sion reduction of special second-order and general higher-order systems via 
dimension reduction of corresponding first-order systems. In Section EJ we 
review the concepts of block-Krylov subspaces and Pade-type reduced-order 
models. In Sectional we introduce the notion of Hermitian higher-order linear 
dynamical systems, and we establish an enhanced Pade-type approximation 
property in the Hermitian case. In Section |H1 we present the SPRIM algo- 
rithm for special second-order systems. In Sectional we report results of some 
numerical experiments with the SPRIM algorithm. Finally, in Section we 
mention some open problems and make some concluding remarks. 

Throughout this paper the following notation is used. Unless stated oth- 
erwise, all vectors and matrices are allowed to have real or complex entries. 
For a complex number a or a complex matrix M, we denote its complex con- 
jugate by a or M, respectively. For a matrix M = [rrijk ], M T := [m/y ] is the 
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transpose of M, and M H := M = [ mkj ] is the conjugate transpose of M. 
For a square matrix P, we write P y if P = P H is Hermitian and if P is 
positive scmidefinite, i.e., x H Px > for all vectors x of suitable dimension. 
We write P y if P = P H is positive definite, i.e., x H Px > for all vectors 
x, except x = 0. The n x n identity matrix is denoted by I n and the zero 
matrix by 0. If the dimension of /„ is apparent from the context, we drop the 
index and simply use /. The actual dimension of will always be clear from 
the context. The sets of real and complex numbers are denoted by M and C, 
respectively. 

2 RCL circuits as first-order and second-order systems 

An important class of electronic circuits is linear RCL circuits that contain 
only resistors, capacitors, and inductors. For example, such RCL circuits are 
used to model the interconnect of VLSI circuits; see, e.g., |H llfil l2*2| . In this 
section, we briefly review the RCL circuit equations and their formulations as 
first-order and second-order linear dynamical systems. 

2.1 RCL circuit equations 

General electronic circuits are usually modeled as networks whose branches 
correspond to the circuit elements and whose nodes correspond to the inter- 
connections of the circuit elements; see, e.g., Such networks are char- 
acterized by Kirchhoff's current law (KCL), Kirchhoff's voltage law (KVL), 
and the branch constitutive relations (BCRs). The Kirchhoff laws depend only 
on the interconnections of the circuit elements, while the BCRs characterize 
the actual elements. For example, the BCR of a linear resistor is Ohm's law. 
The BCRs are linear equations for simple devices, such as linear resistors, 
capacitors, and inductors, and they are nonlinear equations for more complex 
devices, such as diodes and transistors. 

The connectivity of such a network can be captured using a directional 
graph. More precisely, the nodes of the graph correspond to the nodes of the 
circuit, and the edges of the graph correspond to each of the circuit elements. 
An arbitrary direction is assigned to graph edges, so one can distinguish be- 
tween the source and destination nodes. The adjacency matrix, A, of the 
directional graph describes the connectivity of a circuit. Each row of A corre- 
sponds to a graph edge and, therefore, to a circuit element. Each column of 
A corresponds to a graph or circuit node. The column corresponding to the 
datum (ground) node of the circuit is omitted in order to remove redundancy. 
By convention, a row of A contains +1 in the column corresponding to the 
source node, —1 in the column corresponding to the destination node, and 
everywhere else. Kirchhoff's laws can be expressed in terms of A as follows: 

KCL: A T i b = 0, 

(1) 

KVL: Av n =v b . 
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Here, the vectors if, and Vf, contain the branch currents and voltages, respec- 
tively, and v n the non-datum node voltages. 

We now restrict ourselves to linear RCL circuits, and for simplicity, we 
assume that the circuit is excited only by current sources. In this case, A, Vb, 
and ib can be partitioned according to circuit-element types as follows: 



A 



A r 
Ai 



Vb = Vb(t) 



V 9 
Vc 
Vl 



ib = ib{t) 



(2) 



Here, the subscripts i, g, c, and I stand for branches containing current sources, 
resistors, capacitors, and inductors, respectively. Using ©, the KCL and KVL 
equations Q take on the following form: 



AiV ri 



AgV n 



Alii + A T g i g + A T c i c + Ajk = 0, 
v c , Aiv n = Vl. 
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Furthermore, the BCRs can be stated as follows: 



-I(t), 



Gv. 



C dt V - 



T d ■ 



(3) 



(4) 



Here, I{t) is the vector of current-source values, G >~ and C >- are diagonal 
matrices whose diagonal entries are the conductance and capacitance values 
of the resistors and capacitors, respectively, and L y is the inductance 
matrix. In the absence of inductive coupling, L is also a diagonal matrix, but 
in general, L is a full matrix. However, an important special case is inductance 
matrices L whose inverse, the so-called susceptance matrix, S = L^ 1 is sparse; 
see [23 Ell- 
Equations @ and ifUl. together with initial conditions for v n (to) and ii(to) 
at some initial time to, provide a complete description of a given RCL circuit. 
For simplicity, in the following we assume to — with zero initial conditions: 



v n (0) = and ii(0) = 0. 



(5) 



Instead of solving J^J and (Jlj directly, one usually hrst eliminates as many vari- 
ables as possible; this procedure is called modified nodal analysis (MNA) |15l 
I25| . More precisely, using the last three equations in J3J and the first three 
equations in (0J, one can eliminate v g , v c , vi, ii, i g , i c , and is left with the 
coupled equations 



Afl(t) = A T a GA g v n + A T c CA c -v n + Aji u 



(6) 



Aiv n = L—i[ 
at 



for v n and i;. Note that the equations (JSJ) are completed by the initial condi- 
tions (JSJ. 
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For later use, we remark that the energy supplied to the RCL circuit by 
the current sources is given by 



E(t)= / Mr)) T /(T)dT. 



(7) 



Recall that the entries of the vector m are the voltages at the current sources. 
In view of the second equation in (J3J), m is connected to v n by the output 
relation 

Vi = A l v n . (8) 



2.2 RCL circuits as first-order systems 

The RCL circuit equations JfjJ and (JSJ can be viewed as a first-order time- 
invariant linear dynamical system with state vector 



z(t) := 



v n (t) 



and input and output vectors 

u(t):=I(t) and y(t):=v,(t), 
respectively. Indeed, the equations © and JSJ) are equivalent to 

S —z(t) -Az(t) = Bu{t), 
y(t)=B T z(t), 



where 
£ := 







A 



—A^GAg —Af 
Ai 



B 







(9) 



(10) 



(11) 



Note that (|10|l is a system of differential- algebraic equations (DAEs) of first 
order. Furthermore, in view of J^J, the energy Q, which is supplied to the 
RCL circuit by the current sources, is just the integral 



E(t)= f\y(r)) T u(r)dT 
Jo 



(12) 



of the inner product of the input and output vectors of IjlQI) . RCL circuits are 
passive systems, i.e., they do not generate energy, and (|12f> is an important 
formula for the proper treatment of passivity; see, e.g., |21 119j. 
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2.3 RCL circuits as second-order systems 

In this subsection, we assume that the inductance matrix L of the RCL circuit 
is nonsingular. In this case, the RCL circuit equations © and (jSJl can also be 
viewed as a second-order time-invariant linear dynamical system with state 
vector 

x{t) := v n (t), 

and the same input and output vectors © as before. Indeed, by integrating 
the second equation of © and using (J5J, we obtain 

Lii(t) =Aif v n {r)dT. (13) 
Jo 

Since L is assumed to be nonsingular, we can employ the relation l|13(l to 
eliminate ii in J^J. The resulting equation, combined with (JHJ, can be written 
as follows: 

d r* 

P 1 — X {t) + P x(t)+P- 1 / x(T)dT = B U (t), 

dt J (14) 

y(t) = B T x{t). 

Here, we have set 

P X :=A T C CA C , P :— A^GA g , P_ : := Aj A u B:=Aj. (15) 

Note that the first equation in ljl4|) is a system of integro-DAEs. We will 
refer to (|14fl as a special second-order time-invariant linear dynamical system. 
We remark that the input and output vectors of l|14[) are the same as in 
the first-order formulation (|10|) . In particular, the important formula l|12|) for 
the energy supplied to the system remains valid for the special second-order 
formulation H10() . 

If the input vector u(t) is differentiable, then by differentiating the first 
equation of 114|) we obtain the "true" second-order formulation 

p ^ x v +p 4t x v +p -i x v = B i u{t) > (i6) 

y{t) = B T x(t). 

However, besides the additional assumption on the differentiability of u(t), 
the formulation l|ltj|) also has the disadvantage that the energy supplied to 
the system is no longer given by the integral of the inner product of the input 
and output vectors 

u{t) ■= 4 W W and : = Vit) 

at 

of H16|l . Related to this lack of a formula of type H12J) is the fact that the 
transfer function of (|l(jfl is no longer positive real. For these reasons, we prefer 
to use the special second-order formulation (|14|) , rather than the more common 
formulation l|16|) . 
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3 Higher-order linear dynamical systems 

In this section, we discuss our general framework for special second-order and 
higher-order linear dynamical systems. We denote by m and p the number of 
inputs and outputs, respectively, and by I the order of such systems. In the 
following, the only assumption on m, p, and I is that m, p, I > 1. 

3.1 Special second-order systems 

A special second-order m-input p-output time-invariant linear dynamical sys- 
tem of order I is a system of integro-DAEs of the following form: 

d r* 

P 1 —x{t)+P x(t)+P-i / x(r)dT = Bu(t), 
dt J to 

y(t) = Du(t) + Lx(t), ( 17 ) 
x(t ) = x . 

Here, P_i, P , Pi E C NxN , B E C Nxm , D 6 C pxm , and L E C pxN are given 
matrices, t E K is a given initial time, and xq £ C N is a given vector of initial 
values. We assume that the N x N matrix 

sPi+P + -P_i 
s 

is singular only for finitely many values of s E C. 

The frequency-domain transfer function of Ijl7(l is given by 

H{s) = D + l(sP 1 +P + ^P- 1 ^ 1 B. (18) 

Note that 

H :C^(CUoo) pxm 

is a matrix-valued rational function. 

In practical applications, such as the case of RCL circuits described in 
Section[21 the matrices Po and Pi are usually sparse. The matrix P_i, however, 
may be dense, but has a sparse representation of the form 

P_i = FiGFf (19) 

or 

P_i = PiG _1 P 2 H , with nonsingular G, (20) 

where Fx, F2 € C JVxAr ° and G € C JVn><JVo are sparse matrices. We stress that 
in the case <|19[) . the matrix G is not required to be nonsingular. In particular, 
for any matrix P_i G C NxN , there is always the trivial factorization 119|) 
with Pi = F2 — I and G = P_i. Therefore, without loss of generality, in the 
following, we assume that the matrix P_i in l|17fl is given by a product of the 
form (S) or QQ|. 
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3.2 General higher-order systems 

An m-input p-output time-invariant linear dynamical system of order I is a 
system of DAEs of the following form: 

P ^x(t) + Pj_i ^rrx(t) + ■ ■ ■ + P x j t x(t) + P x(t) = Bu(t), 

d 1 - 1 d 
y(t) = D u{t) + L ; _! -TaZix{t) + ■ ■ ■ + L x —x[t) + L x(t). 

Here, P E C NxN , < i < I, B € C Nxm , D G C pxm , and Lj G C pxN , 
< j < I, are given matrices, and N is called the state-space dimension of (|21|l . 
Moreover, in (|21|l . u : [to, oo) C m is a given input function, to G K is a given 
initial time, the components of the vector- valued function x : [to, oo) i— > C w 
are the so-called state variables, and y : [to, oo) i— > C p is the output function. 
The system is completed by initial conditions of the form 



-r-x(t) 

dp y ' 



xf , 0<i<l, (22) 



t=t 



where x ^ G C^, < i < I, are given vectors. 

The frequency-domain transfer function of l|21|) is given by 

H(s) :=D + L(s) (P(s)) ~ X B, s G C, (23) 



where 
and 

Note that 



P(s) := s l Pi + s 1 - 1 ^-! + --- + sP 1 +P (24) 
L(s) := s 1 - 1 ^ + s l - 2 Li- 2 + • ■ ■ + sL x + L . 



P:C^C NxN and L : C i-> C pxN 
are matrix- valued polynomials, and that 

P:C^(CUoo) pxm 

is a matrix- valued rational function. We assume that the polynomial (|24|l . P, 
is regular, that is, the matrix P(s) is singular only for finitely many values of 
s G C; see, e.g., ^] Part II]. This guarantees that the transfer function 
has only finitely many poles. 



3.3 First-order systems 



For the special case 1 = 1, systems of the form l|2H are called first-order 
systems. In the following, we use calligraphic letters for the data matrices and 
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z for the vector of state-space variables of first-order systems. More precisely, 
we always write first-order systems in the form 



S—z(t) -Az(t) =Bu(t), 
at 

y(t) =Vu{t) + Cz(t), 
z(t Q ) = z . 

Note that the transfer function of (|25H is given by 

H(s) = V + £(s£-Ay 1 B. 



(25) 



(26) 



4 Equivalent first-order systems 

A standard approach to treat higher-order systems is to rewrite them as equiv- 
alent first-order systems. In this section, we present such equivalent first-order 
formulations of special second-order and general higher-order systems. 

4.1 The case of special second-order systems 

We start with special second-order systems <|17[) . and we distinguish the two 
cases (Uni and (TZUj) . 

First assume that P-\ is given by l|19|) . In this case, we set 



zx(t):=x{t) and z 2 (t) := F 2 I x(t) dr 



(27) 



to 



By <|19t and (I27L the first relation in (|17(l can be rewritten as follows: 



Pi j f zx{t) + P z 1 (t)+F 1 Gz 2 (t) = Bu(t). 



Moreover, l|27|) implies that 



G H -z 2 (t) - (F 2 G) H Zl (t). 
at 



(28) 



(29) 



It follows from l|27(l - (|29|l that the special second-order system ()17(l (with P_i 
given by H19|) ~) is equivalent to a first-order system l|25|l where 



Zl (t) 

*2(«) 



z(t) 
V:=D, A 



z := 



■>•() 




-Po -F X G 
(F 2 G) H 



C:=[L 0], B:= 
£ : 



Pi 
G H 



(30) 
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The state-space dimension of this first-order system is Ni := N + N , where 
N and N denote the dimensions of Fx G C"^ and G £ C N ° xN °. Note 
that (|26[1 is the corresponding representation of the transfer function IjlH)). H, 
in terms of the data matrices defined in i|3U[l . 

Next, we assume that P_i is given by (|20() . We set 

21 (t) := a;(i) and z 2 (£) := G~ 1 F 2 H [ x(r)d,T. 



In 



The first relation in l|17|) can then be rewritten as 



Pi—z 1 (t)+Poz 1 {t)+F 1 z 2 (t)=Bu(t). 



Moreover, we have 



G-z 2 (t) = F 2 H Zl (t). 



It follows that the special second-order system (|17fl (with P_i given by l|2U|)l 
is equivalent to a first-order system (12311 where 



,{t) := 



X> 



2l (f) 



z := 



x 




£:=[L 0], B:= 



-Po -Pi 
P 2 ff 



£ 



Pi 
G 



(31) 



The state-space dimension of this first-order system is again N% :~ N + No. 
Note that (|26|l is the corresponding representation of the transfer func- 
tion ijl8|l . H, in terms of the data matrices defined in l|31|l . 



4.2 The case of general higher-order systems 

It is well known (see, e.g., ^] Chapter 7]) that any l-th order system with 
state-space dimension N is equivalent to a first-order system with state-space 
dimension TVi :— IN. Indeed, it is easy to verify that the Z-th order system l|21() 
with initial conditions 122[l is equivalent to the first-order system Ij25(l with 
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z{t) 



x(t) 




r r (0) " 




' " 






T (l) 


, B := 






, z := 





\_ dt i-i x(t) _ 




T ('-i) 




B _ 



C:=[L Q L x 








Pi 



A 








Po 



I 



Pi 





p 2 







Pi 



i-1 



(32) 



We remark that l|26l) is the corresponding representation of the Z-order transfer 
function (|23|) . iJ, in terms of the data matrices defined in l|32fl . 



5 Dimension reduction of equivalent first-order systems 

In this section, we discuss some general concepts of dimension reduction of 
special second-order and general higher-order systems via dimension reduction 
of equivalent first-order systems. 

5.1 General reduced-order models 

We start with general first-order systems l|25|) . For simplicity, from now on 
we always assume zero initial conditions, i.e., zq = in l|25l) . We can then 
drop the initial conditions in (|25|l , and consider first-order systems (|25|) of the 
following form: 

e — z(t)-Az(t) =Bu(t), 

dt (33) 
y(t) =Vu{t)+Lz{t). 

Here, A, £ £ C NlxN \ B x E C NlXm , V e C pxm , and C G C pxNl are given 
matrices. Recall that iVi is the state-space dimension of (13311 . We assume that 
the matrix pencil s£ — A is regular, i.e., the matrix s£ — A is singular only 
for finitely many values of s e C. This guarantees that the transfer function 

of iHSJ, 

H(s) :=V + C(s£ -A)~ 1 B, (34) 

exists. 

A reduced- order model of (|33|) is a system of the same form as (|33|l . but 
with smaller state-space dimension. More precisely, a reduced-order model 
of (|33|l with state-space dimension ri\ (< N\) is a system of the form 
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e p(t) -Am =Bu(t), (35) 

y(t) =Vu(t)+ Cz(t), 

where A, £ £ C" lXni , B € C" lX "\ £> 6 C pxm , and £ € C pxni . Again, we 
assume that the matrix pencil s £ — A is regular. The transfer function of 1|35|) 
is then given by 

H(s):=V + C(s£-Ay 1 B. (36) 

Of course, 1|35|1 only provides a framework for model reduction. The real 
problem, namely the choice of suitable matrices A, £, B, C, £>, and sufficiently 
large reduced state-space dimension n\ still remains to be addressed. 



5.2 Reduction via projection 

A simple, yet very powerful (when combined with Krylov-subspace machinery) 
approach for constructing reduced-order models is projection. Let 

VeC N 1 xn 1 ( 3? ) 

be a given matrix, and set 

A:=V H AV, £:=V H £V, B:=V H B C:=CV, V := V. (38) 

Then, provided that the matrix pencil s £ — A is regular, the system 13511 
with matrices given by (|38|l is a reduced-order model of (|33(l with state-space 
dimension n\. 

A more general approach employs two matrices, 

V, W € C NlXni , 

and two-sided projections of the form 

A:=W H AV, £:=W H £V, B := W H B C:=CV, V := V. 

For example, the PVL algorithm |HJEI can De viewed as a two-sided projection 
method, where the columns of the matrices V and W are the first n% right 
and left Lanczos vectors generated by the nonsymmetric Lanczos process |17| . 

All model-reduction techniques discussed in the remainder of this paper 
are based on projections of the type (|38|l . 

Next, we discuss the application of projections (|38|l to first-order sys- 
tems ()33H that arise as equivalent formulations of special second-order and 
higher-oder linear dynamical systems. Recall from Section0]that such equiv- 
alent first-order systems exhibit certain structures. For general matrices (|37|l . 
V, the projected matrices do not preserve these structures. However, as we 
will show now, these structures are preserved for certain types of matrices V. 
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5.3 Preserving special second-order structure 



In this subsection, we consider special second-order systems l|17[) . where P_i is 
either of the form l|19|) or (|20|1 . Recall that the data matrices of the equivalent 
first-order formulations of l|17|) are defined in l|3U[l . respectively (|31|) . 
Let V be any matrix of the block form 



V 



Vi 






v-2 



Vi G c 



iV xra 



(39) 



such that the matrix 



G := V^GV 2 is nonsingular. 



First, consider the case of matrices P_i of the form l|19|l . Using (|3(J[1 and (|39[1 . 
one readily verifies that in this case, the projected matrices l|38[) are as follows: 



A 



-Po 
(F 2 G) H 



C=[L 0], V 
Here, we have set 



-PiG 


D. 



£ 



Pi 






G H 



B 



P Q := VfPoVi, A := V^PiFi, 



B := V^P, 



i := LVi . 



(40) 



(41) 



and 



Pi := (V 1 H F 1 GV 2 ) G-\ F 2 := (Vf F 2 GV 2 ) G~ l . 

Note that the matrices l|40|l are of the same form as the matrices l|3U|l of the 
first-order formulation (1331) of the original special second-order system <|17[) 
(with P_i of the form H19(l L It follows that the matrices (|40|l define a reduced- 
order model in special second-order form, 



Pl^-x(t)+P x(t) 
at 



P-i 



In 



x(t) dr — B u(t), 

y{t) =Du(t) + Lx(t), 



(42) 



where 



P_! := PiGPf . 

We remark that the state-space dimension of (|42fl is n, where n denotes the 
number of columns of the submatrix V\ in (|3l?|l . 

Next, consider the case of matrices P_i of the form (|20|) . Using i|31|) 
and one readily verifies that in this case, the projected matrices <|3^|) 

are as follows: 



A = 
C=\L 



-Po 
F H 



-Pi 




£ = 



Pi 





G 



B = 



(43) 



V = D. 
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Here, Pq, Pi, B, L are the matrices defined in 1(41(1 . and 

Fx := ^Fi^ 2 , F 2 := V 1 H F 2 V 2 . 

Again, the matrices ((43|l are of the same form as the matrices (|31|) of the first- 
order formulation (|3"3T) of the original special second-order system (|T7|) (with 
P-i of the form 1(201) . It follows that the matrices 1(43(1 define a reduced-order 
model in special second-order form ((42(1 . where 

P_i = FiG^Fi 1 . 



5.4 Preserving general higher-order structure 



We now turn to systems (133(1 that are equivalent first-order formulations of 
general Z-th order linear dynamical systems 1)21(1 . More precisely, we assume 
that the matrices in ((33(1 are the ones defined in ((32(1 . 
Let V be any IN x In matrix of the block form 





S n 








... o ■ 









S n 





... o 




v n = 












> S n G 










'■. 






. 







S n . 





^Nxn 



qH q _ f 



(44) 



Although such matrices appear to be very special, they do arise in connection 
with block-Krylov subspaces and lead to Pade-type reduced-order models; 
see Subsection 16.41 below. The block structure ((44(1 implies that the projected 
matrices (I38H are given by 



A 



-I 




B = 



■•■ 

Po Pi 







-/ 



P-2 






-I 

Pl-l 








p 



(45) 



where 
P 




D 

u n i u m 



C=\Ln U 



L 



i-i _ 



v = v, 



0<i<l, B:= S^B, L 3 := LjS n , < j < I. 



It follows that the matrices 1(45(1 define a reduced-order model in Z-th order 
form, 
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Pi jU(t) + Pi-i 4_T*W + ■ ■ • + A jMt) + Po x(t) = Bu(t), 

d 1 - 1 d (46) 

y{t) = D u(t) + -j^x{t) + ■ ■ ■ + Lx -^x(t) + L x(t), 

of the original Z-th order system (|21[l . We remark that the state-space dimen- 
sion of H4t>|> is n, where n denotes the number of columns of the matrix S n 
in (g3J). 



6 Block-Krylov subspaces and Pade-type models 

In this section, we review the concepts of block-Krylov subspaces and Pade- 
type reduced-order models. 

6.1 Pade-type models 

Let so £ C be any point such that the matrix sq £ — A is nonsingular. Recall 
that the matrix pencil s £ — A is assumed to be regular, and so the matrix 
sq £ — A is nonsingular except for finitely many values of s G C. In practice, 
so G C is chosen such that sq £ — A is nonsingular and at the same time, sq 
is in some sense "close" to a problem-specific relevant frequency range in the 
complex Laplace domain. Furthermore, for systems with real matrices A and 
£ one usually selects so G R in order to avoid complex arithmetic. 

We consider first-order systems of the form and their reduced-order 
models of the form (|35fl . By expanding the transfer function 134f) . H, of the 
original system about s$, we obtain 

H(s) = C(s£- Ay 1 B = c(l + (s-so)M) K 

(47) 

= ]T(-l) l £Af 7e(s-s ) 4 , 

where 

M := (s Q £-Ay 1 £ and K := (s £ - A)~ 1 B. (48) 

Provided that the matrix sq £ — A is nonsingular, we can also expand the 
transfer function QHfiJI. H, of the reduced-order model (|35(l about so- This 
gives 

H(s) =£{s£-A)~ l B 

°° (49) 

i=0 

where 

M := (so^-Ay 1 ^ and K := (s £ - A) 
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We call the reduced-order model l|35f) a Pade-type model (with expansion point 
So) of the original system if the Taylor expansions l|47|l and (|49|l agree in 
a number of leading terms, i.e., 

H(s) = H(s)+O((s-s ) q ) (50) 

for some q — q(A 7 £, B, C, V) > 0. 

Recall that the state-space dimension of the reduced-order model (|35|) is 
n\. If for a given ri\, the matrices A, £, B, C, T> in (|35|) are chosen such that 
q = q{n\) in (|5()|l is optimal, i.e., as large as possible, then the reduced-order 
model (1351) is called a Pade model. All the reduced-order models discussed in 
the remainder of this paper are Pade-type models, but they are not optimal 
in the Pade sense. 

The (matrix-valued) coefficients in the expansions l|47|) and l|49(l are often 
referred to as moments. Strictly speaking, the term "moments" should only 
be used when s = 0; in this case, the Taylor coefficients of the Laplace- 
domain transfer function directly correspond to the moments in time domain. 
However, the use of the term "moments" has become common even in the case 
of general sq. Correspondingly, the property l|50[) is now generally referred to 
as "moment matching" . 

We remark that the moment-matching property H50|l is important for the 
following two reasons. First, for large-scale systems, the matrices A and £ are 
usually sparse, and the dominant computational work for moment-matching 
reduction techniques is the computation of a sparse LU factorization of the 
matrix sq£ ~ A. Note that such a factorization is required already even if 
one only wants to evaluate the transfer function H at the point s - Once 
a sparse LU factorization of Sq £ — A has been generated, moments can be 
computed cheaply. Indeed, in view of l|47|) and l|48|) . only sparse back solves, 
sparse matrix products (with £ ), and vector operations are required. Second, 
the moment-matching property l|50[) is inherently connected to block-Krylov 
subspaces. In particular, Pade-type reduced-order models can be computed 
easily be combining Krylov-subspace machinery and projection techniques. In 
the remainder of the section, we describe this connection with block-Krylov 
subspaces. 

6.2 Block-Krylov subspaces 

In this subsection, we review the concept of block-Krylov subspaces induced 
by the matrices M and K defined in Recall that A, £ £ C NlxNl and 

B eC NlXm . Thus we have 

MeC NlxNl and KeC NlXm . (51) 

Next, consider the infinite block-Krylov matrix 



[K Mil M 2 TZ 



(52) 
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In view of l|51|l . the columns of the matrix i|52|) are vectors in C *, and so 
only at most Ni of these vectors are linearly independent. By scanning the 
columns of the matrix l|52fl from left to right and deleting each column that 
is linearly dependent on columns to its left, one obtains the so-called deflated 
finite block-Krylov matrix 

[KM MKM M 2 K^ ••• M**"*- 1 H<*~J } , (53) 

where each block K^> is a subblock of TZ^^ X \ j = 1,2, . . . , j ma x, and K^ :— 
1Z. Let nij denote the number of columns of the j-th block K- ' . Note that 
by construction, the matrix (|53|l has full column rank. The n-th block-Krylov 
subspace (induced by M. and 1Z) JC n (M , 1Z) is defined as the subspace of C^ 1 
spanned by the first n columns of the matrix (|53j) : see, pQ for more details of 
this construction. 

Here, we will only use those block-Krylov subspaces that correspond to 
the end of the blocks in l|53fl . More precisely, let n be of the form 

n = n(j) := mi + + ■ • • + rrij, where 1 < j ' < j max . (54) 

In view of the above construction, the n-th block-Krylov subspace is given by 

K n (M,K) = range [K^ MTZ^ M 2 K^ ••• M^K^]. (55) 

6.3 The projection theorem revisited 

It is well known that the projection approach described in Subsection 15.21 
generates Pade-type reduced-order models, provided that the matrix V in 13 7|) 
is chosen as a basis matrix for the block-Krylov subspaces induced by the 
matrices A4 and 1Z defined in (|48|l . This result is called the projection theorem, 
and it goes back to at least It was also established in (201 1211 122] in 
connection with the PRIMA reduction approach; see JU] for more details. 

One key insight to obtain structure-preserving Pade-type reduced-order 
models via block-Krylov subspaces and projection is the fact that the projec- 
tion theorem remains valid when the above assumption on V is replaced by 
the weaker condition 

K, n (M,K) C range V„. (56) 

In this subsection, we present an extension of the projection theorem (as 
stated in ^U]) to the case (O- 

From now on, we always assume that n is an integer of the form (|54() and 
that 

V n eC NlXni (57) 

is a matrix satisfying l)56[l. Note that (|56|l implies n\ > n. We stress that we 
make no further assumptions about n\. We consider projected models given 
by H38fl with V = V„. In order to indicate the dependence on the dimension 
n of the block-Krylov subspace in 1)56(1. we use the notation 
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A n := V*A V ni £n := V* 5 V„, B„ := V*B, 

(58) 

£„:=£V„, £> n :=X> 

for the matrices defining the projected reduced-order model. In addition 
to (|5fcip . we also assume that the matrix pencil s£ n — An is regular, and 
that at the expansion point sq, the matrix sq £ n — A n is nonsingular. Then 
the reduced-order transfer function 



-1 

T> 

(59) 



£„ (/ + (s - s )M„) 7^ 



= ^(-l) 4 £„^ l ^„( s - So ) 1 

is a well-defined rational function. Here, we have set 

M„ := (s £n — A^) £ n and TZ n :— (s Q £ n - An) B n . (60) 

We remark that the regularity of the matrix pencil s £ n — A n implies that the 
matrix V„ must have full column rank. 

After these preliminaries, the extension of the projection theorem can be 
stated as follows. 

Theorem 1. Let n = n(J) be of the form S5m, and let V n be a matrix satis- 
fying (|56]1 . Then the reduced- order model defined by is a Pade-type model 
with 

Hn( S )=H( S ) + O((s- S0 y). (61) 

Proof. In view of J37JI and the claim i|6T|l holds true if 

M i Tl = V n M n K» for all i = 0, 1, . . . ,j - 1, (62) 

and thus it is sufficient to show H62[) . 

By and l(5^|) . for each i = 0, 1, . . . , j — 1, there exists a matrix pi such 
that 

M l n = V n p i . (63) 

Moreover, since V n has full column rank, each matrix pi is unique. In fact, we 
will show that the matrices pi in (|63() are given by 

Pi =M n K n , , 0.1 ./ 1. (64) 

The claim Il(i2l) then follows by inserting l|64|l into l|63|) . 

We prove l|64|l by induction on i. Let i = 0. In view of (|48|l and i|63[l . we 

have 

VnPo =n= (s a £-Ay 1 B. (65) 
Multiplying l|65|) from the left by 
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(a £ n -A n )~ 1 V^(8o£ -A) (66) 

and using the definition of 1Z n in (|60|l . it follows that po = TZ n - This is just 
the relation l|64|) for i = 0. 

Now let 1 < i < j — 1, and assume that 

Pi-i =M l n - 1 ^„. (67) 

Recall that pi-i satisfies the equation (|63(l (with i replaced by i — 1), and thus 
we have .A/P -1 7?. = V„ Together with (ESI and (jHZI), it follows that 

V nPl = M l K = M{M l - 1 K) =M{V nPi - X ) = MV^M); 1 H n ) . (68) 

Note that, in view of the definition of M. in (|48() . we have 

V* (s £ - A)M V n = V%£ V n = S n . (69) 

Multiplying (|68(l from the left by the matrix (j66(l and using <69|l as well as 
the definition of M. n in 1)60(1 . we obtain 

Pi = (s £n - AO "^n (M'n 1 TZ n ) = M\ K n . 

Thus the proof is complete. 



6.4 Structure-preserving Pade-type models 

We now turn to structure-preserving Pade-type models. Recall that, in Sub- 
sections 15.31 and 15.41 we have shown how special second-order and general 
higher-order structure is preserved by choosing projection matrices of the 
form H39J1 and (|44J) . respectively. Moreover, in Subsection 16 . ?>\ we pointed out 
that projected models are Pade-type models if (|5l)|l is satisfied. It follows 
that the reduced-order models given by the projected data matrices lf5"%l) are 
structure-preserving Pade-type models, provided that the matrix V„ in 1571) 
is of the form l|39|l . respectively l|44(l . and at the same time fulfills the condi- 
tion 1)56(1. Next we show how to construct such matrices V n - 
Let 

V n eC NlXn (70) 

be any matrix whose columns span the rt-th block-Krylov subspace K. n , 1Z) , 
i.e., 

K. n (M,n) = range V„. (71) 

First, consider the case of special second-order systems l)17|l . where P-\ is 
either of the form (|19|) or 120fl . In this case, we partition V n as follows: 

Fi e c Nxn , v 2 e c NoXn . (72) 



V n = 



Vi 
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Using the blocks in (|72|l . we set 



V„ := 



Vi 

v 2 



(73) 



Clearly, the matrix (|73|) is of the form (|39[) . and thus the projected mod- 
els generated with V n preserve the special second-order structure. Moreover, 
from JHJ-O, it follows that 

fC n (M,Tl) = range V„ C range V n , 

and so condition (|56|l is satisfied. Thus, the projected models are Pade-type 
models and preserve second-order structure. 

Next, we turn to the case of general higher-order systems (|21|l . In 12 , 
we have shown that in this case, the block-Krylov subspaces induced by the 
matrices M. and 1Z, which are given by l|32|) and (|48|l . exhibit a very special 
structure. More precisely, the n-dimensional subspace K, n (M,lZ) of C lN can 
be viewed as I copies of an n-dimensional subspace of C^. Let S n 6 C Nxn 
be a matrix whose columns form an orthonormal basis of this n-dimensional 
subspace of C N , and set 





S n 










Sn 


V„ := 










. 














(74) 



From the above structure of the n-dimensional subspace JC n {M. , 1Z) , it fol- 
lows that V„ satisfies the condition (|56|) . Furthermore, the matrix V„ is of 
the form (|44(l . Thus, the projected models generated with V n are Pade-type 
models and preserve higher-order structure. 

In the remainder of this paper, we assume that V n are matrices given 
by l|73|l in the case of special second-order systems, respectively Ij74(l in the 
case of higher-order systems, and we consider the corresponding structure- 
preserving reduced-order models with data matrices given by (|58|l . 



7 Higher accuracy in the Hermitian case 

For the structure-preserving Pade-type models introduced in Subsection 16.41 
the result of Theorem ^ can be improved further, provided the underlying 
special second-order or higher-order linear dynamical system is Hermitian, 
and the expansion point sq is real, i.e., 



s e M. 



(75) 
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More precisely, in the Hermitian case, the Pade-type models obtained via V„ 
match 2j(n) moments, instead of just j(n) in the general case; see Theorem^ 
below. We remark that for the special case of real symmetric second-order 
systems and expansion point s = 0, this result can be traced back to [23] ■ 

In this section, we first give an exact definition of Hermitian special second- 
order systems and higher-order systems, and then we prove the stronger 
moment- matching property stated in Theorem 



7.1 Hermitian special second-order systems 



We say that a special second-order system (|17|) is Hermitian if the matrices 
in (|17fl and 119|) . respectively (|20|l . satisfy the following properties: 



L = B 



H 



Pn = P 



H 



P 1 =Pl I , Fi=F 2 , G = G 



H 



(76) 



Recall that RCL circuits are described by special second-order systems of 
the form 1)14(1 with real matrices defined in l(15JI . Clearly, these systems are 
Hermitian. 

Using H75fl . (|76fl . and (|19fl . respectively (|2U[) . one readily verifies that the 
data matrices (|30|l . respectively (|31|l . of the equivalent first-order formula- 
tion ll.').''>[l satisfy the relations 



J(s £-A) = (s Q £ -A) H J, J£ = £J, J = J H , 
C H = JB, 



(77) 



where 



J 



In 
-I No 



Since the reduced-order model is structure-preserving, the data matrices (|58|) 
satisfy analogous relations. More precisely, we have 



C H — 7 E 



(78) 



where 



In 






-In 



7.2 Hermitian higher-order systems 

We say that a higher-order system l|21[l is Hermitian if the matrices in H21|) 
satisfy the following properties: 

Pt = P", 0<i<l, L = B H , Lj = 0, 1 < j < Z — 1. (79) 
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In this case, we define matrices 



1-3 
i=0 



+») 



and set 



J:= 



-s I 



I 

•■• 




sol 



/ 





s I 



j = 0, 



Pi P2 



Pl-l ■■' 

Pi o 



I 

Note that, in view of we have 

P 3 =Pf, J =0,1,...,/ 



3 
: 





(80) 



(81) 



Using (J75J|-|jHlJ, one can verify that the data matrices A, £, B, C defined 
in (|32[1 satisfy the following relations: 

J(s £-A) = (s £-A) H J, J£ = £ H J, C H = JB. (82) 

Since the reduced-order model is structure-preserving, the data matrices Q58JI 
satisfy the same relations. More precisely, we have 



J7n (^0 £ n An^) \Sq £ n An) 3m *Jn £n — £n ^n: 

C H — 7 B 
where J n is defined in analogy to J . 



(83) 



7.3 Key relations 

Our proof of the enhanced moment-matching property in the Hermitian case 
is based on some key relations that hold true for both special second-order 
and higher-order systems. In this subsection, we state these key relations. 

Recall the definition of the matrix M. in l|48|) . The relations J77J), respec- 
tively l|82J) . readily imply the following identity: 



M H J = J£(s Q £ -A) \ 

It follows from (H2J that 

(M H Y J = j(e (s £ -A)' 1 )', i = 0,1, 
Similarly, the relations f78"| . respectively (JEHJ), imply 



(84) 



(85) 
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It follows that 

{M*Y J = J n (s n (s S n -A n y i y, t = 0, 1 (86) 

Also, recall from J77J, respectively (|S2|l . that 

C H = JB, (87) 
and from (JTSJ), respectively (JSSJ), that 

£%=J n B n . (88) 
Finally, one readily verifies the following relation: 

V ? f J£ V„ =Jn£n. (89) 

7.4 Matching twice as many moments 

In this subsection, we present our enhanced version of Theorem^for the case 
of Hcrmitian special second-order or higher-order systems. 
First, we establish the following proposition. 

Proposition 1. Let n = n(j) be of the form ]54[ - Then, the data matri- 
ces <|5#f) of the structure-preserving Pade-type model satisfy 

CM l V n = C n AV n for all i = 0, 1, . . . , j. (90) 

Proof. Recall that C n = C V„. Thus holds true for i = 0. 
Let 1 < i < j. In view of l|85[) . we have 

(M H yj = j(£{s s-Ay i y. 

Together with l|g7J). it follows that 

{M H YC H ={M H YJB = J (S ( 8Q 8- AY^B. 

Since (so £ -A) B = K, it follows that 

(M H y C H = J£ ((s £- A)~ l £y~ 1 K = J £ M^K. 
Using (with i replaced by i - 1), JH§J|, JSSJ, and (jSHJ, we obtain 
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Vj? (M H ) l C H = V*J£ (M^K) 
= V% J£V n M^ 1 K n 
= (V* J£Vn)(Mt 1 Kn) 

= J n £ n M i ~ x K n 

= J n £ n M)- 1 (so£-AY 1 B n 

= J n (£n(s Q £-A)~ i yB n 
= (M»yj n B n =(M*y£%. 

Thus the proof is complete. 

The following theorem contains the main result of this section. 

Theorem 2. Let n = n(j) be of the form \54\ - In the Hermitian case, the 
structure-preserving Pade-type model defined by the data matrices satis- 
fies: 

H n (s) = H(s) + O((s-s ) 2j( - n) ). 

Proof. Let j = j(n). We need to show that 

CM i K= CnMiKn forall i = 0, 1, . . . , 2j - 1. (91) 

By JHlJl and (jHOJ, we have 

£M ll+l2 K= (CM ll )(M Z2 K) 

= (CM^)(V n M\?K n ) 
= (CM^V n )(M\ 2 Tl n ) 
= (C n MH) Kn) = £n MH+^ TZ n 

for all i\ = 0, 1, . . . , j — 1 and i 2 = 0, 1, . . . , j. This is just the desired rela- 
tion (|91J) . and thus the proof is complete. 



8 The SPRIM algorithm 

In this section, we apply the machinery of structure-preserving Pade-type 
reduced-order modeling to the class of Hermitian special second-order systems 
that describe RCL circuits. 

Recall from Section [21 that a first-order formulation of RCL circuit equa- 
tions is given by (|10|) with data matrices defined in (|llf). Here, we consider 
first-order systems (|10|l with data matrices of the slightly more general form 



A = 



-Po -F 
F H 



£ = 



Pi 
G 



B = 



(92) 
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Here, it is assumed that the subblocks Pq, Pi, an d B have the same number 
of rows, and that the subblocks of A and £ satisfy Pq >z 0, P\ >z 0, and 
G >- 0. Note that systems (|1(J|) with matrices ()92[) are in particular Hermitian. 
Furthermore, the transfer function of such systems is given by 

H(s) =B H (s£-Ay 1 B. 

The PRIMA algorithm is a reduction technique for first-order 

systems H10|) with matrices of the form (1921) . PRIMA is a projection method 
that uses suitable basis matrices for the block-Krylov subspaces K. n (A4 , 7Z) ; 
see [5]. More precisely, let V„ be any matrix satisfying (|70|l and (|71|) . The 
corresponding n-th PRIMA model is then given by the projected data matrices 

An := V* A V„, 4 := V*£ V n , B n := V%B. 

The associated transfer function is 

H n ( S ) ^B^sL-AnY 1 ^. 

For n of the form l|54|l . the PRIMA transfer function satisfies 

H(s) = H(s) + O ({s - s ) j(n) ) • (93) 

Recently, we introduced the SPRIM algorithm as a structure-preserving 
and more accurate version of PRIMA. SPRIM employs the matrix V„ obtained 
from V„ via the construction (|72|l and l|73|) . The corresponding n-th SPRIM 
model is then given by the projected data matrices 

A„ := V* AV n , £„ := V*£ V n , B n := V* B. 

The associated transfer function is 

H n ( S ) ^B^S^-Ar^Bn. 

In view of Theorem [21 we have 

H( S )=H(s) + o((s- So ) 2 ^), 

which suggests that SPRIM is "twice" as accurate as PRIMA. 
An outline of the SPRIM algorithm is as follows. 

Algorithm 1 (SPRIM algorithm for special second-order systems) 



Input: matrices 
A = 



-Po 



-F 




£ = 



Pi 
G 



B = 



where the subblocks Pq, Pi, and B have the same number of rows, and the 
subblocks of A and £ satisfy Pq >z 0, Pi >z 0, and G >~ 0; 
an expansion point sq G M. 
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Formally set 



M = (s £ -A)~ X C, K = (s £-A)~ l B. 



• Until n is large enough, run your favorite block-Krylov subspace method 
(applied to M. and TV) to construct the columns of the basis matrix 

V n = [vi v 2 ■■■ v n ] 

of the n-th block-Krylov subspace K. n (A4,lZ^ , i.e., 

span V n = K n (M , TZ) . 

• Let 

Vn 



Vi 



be the partitioning ofV n corresponding to the block sizes of A and £. 
Set 

P Q =V 1 H P 1 V 1 , F = V 1 H FV 2) P x = VFP x Vi, G = V 2 H GV 2 , 
B = VfB, 



-P —F 
F H 



£n = 



Pi 
G 



Output: the reduced-order model H n in first- order form 

H n (s)=B%(s£ n -A n y 1 B n 
and in second-order form 



Hjs) = B H [sP 1 +P + -FG- X F H 
■ s 



B. 



(94) 



(95) 



We remark that the main computational costs of the SPRIM algorithm is 
running the block Krylov subspace method to obtain V n . This is the same as 
for PRIMA. Thus generating the PRIMA reduced-order model H n and the 
SPRIM reduced-order model H n involves the same computational costs. 

On the other hand, when written in first-order form l|94|) . it would ap- 
pear that the SPRIM model has state-space dimension 2n, and thus it would 
be twice as large as the corresponding PRIMA model. However, unlike the 
PRIMA model, the SPRIM model can always be represented in special 
second-order form l)95[l: see Subsection 15.31 In Ij95(l . the matrices Pi, Pq, and 
P_i := FG^F 11 are all of size n x n, and the matrix B is of size n x m. 
These are the same dimensions as in the PRIMA model (I93f) . Therefore, the 
SPRIM model H n (written in second-order form l|95() ) and of the correspond- 
ing PRIMA model H n indeed have the same state-space dimension n. 
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9 Numerical examples 

In this section, we present results of some numerical experiments with the 
SPRIM algorithm for special second-order systems. These results illustrate the 
higher accuracy of the SPRIM reduced-order models vs. the PRIMA reduced- 
order models. 



9.1 A PEEC circuit 

The first example is a circuit resulting from the so-called PEEC discretiza- 
tion of an electromagnetic problem. The circuit is an RCL network con- 
sisting of 2100 capacitors, 172 inductors, 6990 inductive couplings, and a single 
resistive source that drives the circuit. The circuit is formulated as a 2-port. 
We compare the PRIMA and SPRIM models corresponding to the same di- 
mension n of the underlying block Krylov subspace. The expansion point 
s = 2tt x 10 9 was used. In Figure ^ we plot the absolute value of the (2, 1) 
component of the 2 x 2-matrix-valued transfer function over the frequency 
range of interest. The dimension n = 120 was sufficient for SPRIM to match 
the exact transfer function. The corresponding PRIMA model of the same 
dimension, however, has not yet converged to the exact transfer function in 
large parts of the frequency range of interest. Figure ^ clearly illustrates the 
better approximation properties of SPRIM due to matching of twice as many 
moments as PRIMA. 



- 1 Exact 
- - PRIMA model 




Fig. 1. \H 2 ,i | for PEEC circuit 
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9.2 A package model 

The second example is a 64-pin package model used for an RF integrated cir- 
cuit. Only eight of the package pins carry signals, the rest being either unused 
or carrying supply voltages. The package is characterized as a 16-port com- 
ponent (8 exterior and 8 interior terminals). The package model is described 
by approximately 4000 circuit elements, resistors, capacitors, inductors, and 
inductive couplings. We again compare the PRIMA and SPRIM models cor- 
responding to the same dimension n of the underlying block Krylov subspace. 
The expansion point s = 5ir x 10 9 was used. In Figure|21 we plot the absolute 
value of one of the components of the 16 x 16-matrix-valued transfer function 
over the frequency range of interest. The state-space dimension n = 80 was 
sufficient for SPRIM to match the exact transfer function. The corresponding 
PRIMA model of the same dimension, however, does not match the exact 
transfer function very well near the high frequencies; see Figure 




10 

10° 10° 

Frequency (Hz) 



Fig. 2. The package model 



9.3 A mechanical system 

Exploiting the equivalence (see, e.g., J^I) between RCL circuits and mechan- 
ical systems, both PRIMA and SPRIM can also be applied to reduced-order 
modeling of mechanical systems. Such systems arise for example in the mod- 
eling and simulation of MEMS devices. In Figure 01 we show a comparison 
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of PRIMA and SPRIM for a finite-element model of a shaft. The expansion 
point sq = w X 10 3 was used. The dimension n — 15 was sufficient for SPRIM 
to match the exact transfer function in the frequency range of interest. The 
corresponding PRIMA model of the same dimension, however, has not con- 
verged to the exact transfer function in large parts of the frequency range 
of interest. Figure 01 again illustrates the better approximation properties of 
SPRIM due to the matching of twice as many moments as PRIMA. 



10 Concluding remarks 

We have presented a framework for constructing structure-preserving Pade- 
type reduced-order models of higher-order linear dynamical systems. The 
approach employs projection techniques and Krylov-subspace machinery for 
equivalent first-order formulations of the higher-order systems. We have shown 
that in the important case of Hermitian higher-order systems, our structure- 
preserving Pade-type model reduction is twice as accurate as in the general 
case. Despite this higher accuracy, the models produced by our approach are 
still not optimal in the Pade sense. This can be seen easily by comparing 
the degrees of freedom of general higher-order reduced models of prescribed 
state-space dimension, with the number of moments matched by the Pade- 
type models generated by our approach. Therefore, structure-preserving true 
Pade model reduction remains an open problem. 
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Fig. 4. A mechanical system 



Our approach generates reduced models in higher-order form via equiv- 
alent first-order formulations. It would be desirable to have algorithms that 
construct the same reduced-order models in a more direct fashion, without 
the detour via first-order formulations. Another open problem is the "opti- 
mal" way of constructing basis vectors for the structured Krylov subspaces 
that arise for the equivalent first-order formulations. In particular, an algo- 
rithm for this task should be both computationally efficient and numerically 
stable. Some related work on this problem is described in the recent report 
|18j . but many questions remain open. Finally, the proposed approach is a 
projection technique, and as such, it requires the storage of all the vectors 
used in the projection. This clearly becomes an issue for systems with very 
large state-space dimension. 
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